Zadanie 1

#wylosowanie węzłów interpolacji
xs =  1:1:10
ys = [rand() for x in xs]

10-element Array{Float64,1}:
 0.9864235340749192
 0.6575389408765071
 0.6457039436286105
 0.8017304324871597
 0.442640668548663
 0.6568724948239182
 0.9097921885654445
 0.9749497542537617
 0.292194472997251
 0.8862600110142869

using Plots
scatter(xs,ys, label="węzły interpolacji")

zgodnie z wikipedią, wielomian interpolacyjny Lagrange'a ma postać $$ w(x)= \sum_{i=0}^{n}y_i*\prod_{j=0\wedge j\neq i}^{n}\frac{x-x_j}{x_i-x_j} $$

function lagrange_polynomial(point,x,y)
    sum = 0
    for i = 1:length(y)
        prod = 1
        for j = 1:length(x)
            if i != j
                prod = prod * ((point - x[j])/(x[i] - x[j]))
            end
        end
        sum = sum + y[i]*prod
    end
    return sum
end

lagrange = function(point)
    lagrange_polynomial(point,xs,ys)
end

plot!(lagrange, extrema(xs)...,label="interpolacja Lagrange'a", title="przykład interpolacji wielomianem Lagrange'a", xlabel="x",ylabel="f(x)")

Zadanie 2

function generate_diff(x,y)
    #zdefiniowanie macierzy ilorazów różnicowych
    newt = zeros(Float32, length(x), length(x))

    #obliczenie brzegowych wartości macierzy ilorazów:
    for i = 1:length(y)
        newt[i, 1] = y[i]
    end

    #obliczenie pozostałych elementów macierzy

    for j = 2:length(y)
        for i = 1:(length(y)-j+1)
            newt[i, j] = (newt[i+1,j-1]-newt[i, j-1])/(x[j+i-1] - x[i])
        end
    end
    newt
end



function horner(point, i, newt,x,y)
    if i < length(x)
        return newt[1,i]+(point-x[i])*horner(point, i+1, newt,x,y)
    else
        return newt[1, i]
    end
end

newton = generate_diff(xs,ys)

function f(point)
    i = 1
    horner(point, i, newton,xs,ys)
end
plot!(f, extrema(xs)..., title="przykład interpolacji wielomianem Lagrange'a i Newtona", label="interpolacja Newtona")

Zadanie 3

#using Polynomials

from_package=fit(xs, ys)

plot!(from_package, extrema(xs)..., title="przykład interpolacji Lagrange'a, Newtona oraz z pakietu 'Polynomials'", label="interpolacja z pakietu 'Polynomials'")

Można zauważyć, że wszystkie wykresy idealnie nakładają się na siebie. Jest to spowodowane tym, że nie istnieją dwa różne wielomiany stopnia n, które przyjmują takie same wartości w tych samych n+1 punktach

Zadanie 4

using Random
using DataFrames

sleep(0.001)
function test_time()
    times = zeros(Float64, 10000, 6)
    for i = 1:1:1000
        for j = 1:10
            x =  1:1:i
            y = [rand() for element in x]

            times[(i-1)*10+j, 1] = i
            times[(i-1)*10+j, 2] = @elapsed lagrange_polynomial(rand(),x,y)
            times[(i-1)*10+j, 3] = @elapsed begin
                newton1 = generate_diff(x,y)
            end
            times[(i-1)*10+j, 4] = @elapsed horner(rand(),1,newton1,x,y)

            times[(i-1)*10+j, 5] = @elapsed begin
                func = fit(x,y)
            end
            times[(i-1)*10+j, 6] = @elapsed func(rand())
        end

    end
    times
end
result = test_time()

10000×6 Array{Float64,2}:
    1.0  3.01e-7    4.99e-7    3.01e-7    8.0e-6     7.0e-7
    1.0  0.0        9.9e-8     1.0e-7     7.01e-7    1.0e-7
    1.0  9.9e-8     1.0e-7     0.0        3.01e-7    1.0e-7
    1.0  1.0e-7     1.0e-7     0.0        3.0e-7     1.0e-7
    1.0  1.0e-7     0.0        0.0        3.99e-7    1.0e-7
    1.0  1.0e-7     1.01e-7    0.0        4.0e-7     0.0
    1.0  0.0        0.0        0.0        4.01e-7    1.0e-7
    1.0  1.0e-7     9.9e-8     0.0        4.01e-7    0.0
    1.0  0.0        0.0        0.0        3.99e-7    0.0
    1.0  0.0        0.0        0.0        7.99e-7    1.01e-7
    2.0  2.0e-7     3.99e-7    2.01e-7    6.99e-7    1.01e-7
    2.0  0.0        1.0e-7     0.0        3.99e-7    1.0e-7
    2.0  0.0        1.01e-7    1.0e-7     3.99e-7    0.0
    ⋮                                                ⋮
  999.0  0.0031006  0.0022516  1.81e-5    0.0014502  3.0e-7
  999.0  0.0031038  0.0022637  1.8501e-5  0.0014526  3.0e-7
 1000.0  0.003108   0.0025184  1.7801e-5  0.0014561  2.0e-7
 1000.0  0.0031109  0.0023409  1.7e-5     0.0014534  2.0e-7
 1000.0  0.0031278  0.0023133  1.74e-5    0.0014528  3.99e-7
 1000.0  0.003109   0.0022509  1.62e-5    0.0014172  2.01e-7
 1000.0  0.0031045  0.0080851  1.6901e-5  0.0014222  4.01e-7
 1000.0  0.0030303  0.0021998  1.65e-5    0.0014646  2.99e-7
 1000.0  0.0030294  0.0022039  1.6099e-5  0.0014173  2.0e-7
 1000.0  0.0030287  0.0022583  1.5199e-5  0.0014146  2.0e-7
 1000.0  0.0030988  0.0025664  1.5701e-5  0.0014177  3.01e-7
 1000.0  0.0030339  0.0022227  1.53e-5    0.0014288  2.0e-7

df = DataFrame()
df.ilość_punktów = result[:,1]
df.metoda_lagrangea = result[:,2]
df.generowanie_ilorazów_różnicowych = result[:,3]
df.wartości_wielomianu_newtona = result[:,4]
df.tworzenie_wielomianu_z_Polynomials = result[:,5]
df.wartości_z_Polynomials = result[:,6]
df

10,000 rows × 6 columns (omitted printing of 2 columns)

	ilość_punktów	metoda_lagrangea	generowanie_ilorazów_różnicowych	wartości_wielomianu_newtona
	Float64	Float64	Float64	Float64
1	1.0	3.01e-7	4.99e-7	3.01e-7
2	1.0	0.0	9.9e-8	1.0e-7
3	1.0	9.9e-8	1.0e-7	0.0
4	1.0	1.0e-7	1.0e-7	0.0
5	1.0	1.0e-7	0.0	0.0
6	1.0	1.0e-7	1.01e-7	0.0
7	1.0	0.0	0.0	0.0
8	1.0	1.0e-7	9.9e-8	0.0
9	1.0	0.0	0.0	0.0
10	1.0	0.0	0.0	0.0
11	2.0	2.0e-7	3.99e-7	2.01e-7
12	2.0	0.0	1.0e-7	0.0
13	2.0	0.0	1.01e-7	1.0e-7
14	2.0	1.0e-7	1.0e-7	0.0
15	2.0	9.9e-8	2.0e-7	0.0
16	2.0	0.0	1.0e-7	9.9e-8
17	2.0	1.0e-7	1.0e-7	9.9e-8
18	2.0	1.0e-7	1.0e-7	0.0
19	2.0	0.0	1.01e-7	1.0e-7
20	2.0	0.0	0.0	0.0
21	3.0	2.01e-7	2.0e-7	1.0e-7
22	3.0	1.0e-7	1.0e-7	9.9e-8
23	3.0	1.399e-6	1.01e-7	0.0
24	3.0	2.0e-7	1.01e-7	1.0e-7
25	3.0	1.0e-7	1.0e-7	1.01e-7
26	3.0	0.0	1.0e-7	0.0
27	3.0	9.9e-8	1.0e-7	1.0e-7
28	3.0	1.0e-7	9.9e-8	1.0e-7
29	3.0	1.0e-7	9.9e-8	1.0e-7
30	3.0	1.0e-7	0.0	0.0
⋮	⋮	⋮	⋮	⋮

using Statistics
gd =  groupby(df, :ilość_punktów)
avg_sd_df = combine(gd, :metoda_lagrangea => mean,
    :generowanie_ilorazów_różnicowych => mean,
    :wartości_wielomianu_newtona => mean,
    :tworzenie_wielomianu_z_Polynomials => mean, 
    :wartości_z_Polynomials => mean,
    :metoda_lagrangea => std,
    :generowanie_ilorazów_różnicowych => std,
    :wartości_wielomianu_newtona => std,
    :tworzenie_wielomianu_z_Polynomials => std, 
    :wartości_z_Polynomials => std)
avg_sd_df

1,000 rows × 11 columns (omitted printing of 7 columns)

	ilość_punktów	metoda_lagrangea_mean	generowanie_ilorazów_różnicowych_mean	wartości_wielomianu_newtona_mean
	Float64	Float64	Float64	Float64
1	1.0	8.0e-8	9.98e-8	4.01e-8
2	2.0	5.99e-8	1.301e-7	5.99e-8
3	3.0	2.399e-7	1.0e-7	7.0e-8
4	4.0	1.103e-7	1.404e-7	1.001e-7
5	5.0	1.703e-7	1.2e-7	5.99e-8
6	6.0	1.8e-7	1.6e-7	1.0e-7
7	7.0	2.2e-7	2.0e-7	1.1e-7
8	8.0	2.798e-7	2.002e-7	1.202e-7
9	9.0	3.298e-7	2.001e-7	1.502e-7
10	10.0	3.906e-7	2.595e-7	1.5e-7
11	11.0	4.304e-7	3.299e-7	1.798e-7
12	12.0	5.802e-7	3.295e-7	1.703e-7
13	13.0	5.9e-7	3.701e-7	1.999e-7
14	14.0	7.1e-7	4.2e-7	2.1e-7
15	15.0	7.799e-7	3.899e-7	2.0e-7
16	16.0	8.899e-7	4.602e-7	2.101e-7
17	17.0	9.799e-7	5.401e-7	2.499e-7
18	18.0	1.0801e-6	6.098e-7	2.503e-7
19	19.0	1.2101e-6	6.499e-7	2.302e-7
20	20.0	1.3198e-6	7.202e-7	2.997e-7
21	21.0	1.4299e-6	7.596e-7	3.004e-7
22	22.0	1.5202e-6	7.796e-7	3.202e-7
23	23.0	1.6803e-6	1.07e-6	3.196e-7
24	24.0	1.8299e-6	1.01e-6	3.306e-7
25	25.0	2.0299e-6	1.3599e-6	4.397e-7
26	26.0	2.2199e-6	1.37e-6	3.402e-7
27	27.0	2.3196e-6	1.3002e-6	3.301e-7
28	28.0	2.4701e-6	1.4697e-6	4.002e-7
29	29.0	2.66e-6	1.4903e-6	3.897e-7
30	30.0	2.85e-6	1.6604e-6	4.0e-7
⋮	⋮	⋮	⋮	⋮

scatter([avg_sd_df.metoda_lagrangea_mean avg_sd_df.generowanie_ilorazów_różnicowych_mean avg_sd_df.wartości_wielomianu_newtona_mean avg_sd_df.tworzenie_wielomianu_z_Polynomials_mean avg_sd_df.wartości_z_Polynomials_mean], layout = 5, title=["Lagrange'a" "ilorazy różnicowe" "Newtona" "Polynomials 1" "Polynomials 2"],
        xlabel="ilość punktów", ylabel="czas[s]", yerror=[avg_sd_df.metoda_lagrangea_std 
        avg_sd_df.generowanie_ilorazów_różnicowych_std
        avg_sd_df.wartości_wielomianu_newtona_std
        avg_sd_df.tworzenie_wielomianu_z_Polynomials_std
        avg_sd_df.wartości_z_Polynomials_std])

niestety z powodu małej ilości miejsca, musiałem skrócić nazwy wykresów. Polynomial 1 oznacza wykres, na którym podane są czasy tworzenia wielomianu za pomocą pakietu "Polynomials", a Polynomials 2 oznacza czasy obliczania wielomianu utworzonego wcześniej z pakietu "Polynomials"

Zadanie 5

function neville(x,y,point)
    result = zeros(Float32, length(x))
    for k = 1:(length(x))
        result[k] = y[k];   
        for j in k-1:-1:1
            result[j] = result[j+1]+(result[j + 1] - result[j])*(point - x[k]) / (x[k] - x[j]);
        end
    end
    result[1]
end

function test_neville(x,y)
    for i = 1:length(x)
        val = rand()*10
        print(lagrange_polynomial(val,x,y), " ",neville(x,y,val),"\n")
    end
end

test_neville(xs,ys)

-6.875124267197301 -6.8751254
0.9424087926772629 0.9424089
0.5047162419831147 0.5047163
0.4606721214728949 0.4606721
0.4708525599649966 0.47085255
0.28948811404637315 0.28948903
0.8771440217058487 0.87714416
0.8418664214695775 0.84186643
0.7587729660259565 0.75877297
0.35513185855755935 0.35513195

wartości wielomianu w algorytmie Nevilla i Lagrange'a są takie same

Zadanie 6

#using Pkg
#Pkg.add("Interpolations")
#using 

xsf=1:0.01:10
itp = interpolate(ys, BSpline(Cubic(Free(OnGrid()))))
ysf3=[itp(x) for x in xsf]
plot!(xsf,ysf3, label="BSpline cubic interpolation grid free", title="różne przykłady interpolacji")
#itp = interpolate(ys, BSpline(Cubic(Flat(OnGrid()))))

interp_quad = interpolate(ys, BSpline(Quadratic(Natural(OnGrid()))))
ysf3=[itp(x) for x in xsf]
plot!(xsf,ysf3, label="BSpline quadratic interpolation on Natural grid")

Widać, że BSpline'y (zarówno sześcienne jak i kwadratowe) są pozbawione efektu Runghego na końcach przedziałów. Efekt ten występuje przy interpolacji wielomianowej przy końcach przedziałów i znacząco obniża dokładność przybliżeń funkcji interpolacją wielomianową. Jak widać, jednym ze sposobów na pozbycie się tego efektu jest użycie funkcji sklejanych (np. BSpline'ów)